#################################
#################################
#################################
#
# Replication code for:
#
# "Until the Bitter End? The Diffusion of Surrender Across Battles"
#
# International Organization (forthcoming)
# This version: 9 October, 2017
#
# Authors:
# Todd Lehmann (tlehmann@umich.edu)
# Department of Political Science
# University of Michigan
#
# Yuri Zhukov (zhukov@umich.edu, corresponding)
# Department of Political Science
# University of Michigan
#
#################################
#################################
#################################


# Set working directory
setwd("~/Downloads/ReplicationArchive/")

#################################
## Set up data
#################################

# Load functions and data
rm(list=ls())
source("functions.R")
load("Data/data_IO.RData")

# Ground warfare subset
dataX0 <- dataX[dataX$SEA_BATTLE==0&dataX$AIR_BATTLE==0,]
locations.poly0 <- locations.poly[which(locations.poly$BATTLE_ID%in%dataX0$BATTLE_ID),]

# dependent/independent variables
dvarz1 <- c("ldr_surrendered","lpow_max","lkia_max","ltotal_cas1","ler1","ller1")
xvarz1 <- c("dplus","cinc","opp_cinc","diff_cinc","xpolity","opp_xpolity","diff_xpolity","polity2","opp_polity2","diff_polity2","lpow_max","lkia_max","ltotal_cas1","ldeploy_dist","diff_polity2_dum","diff_polity2_dum_alt","gr_polity2","gr_polity2_new","gr_polity2_alt","gr_xpolity","gr_xpolity_new","gr_xconst","battle_size","size_ratio","lsize_ratio_1","lsize_ratio_2","lsize_ratio_3"
)


################################
################################
## Table 1: summary statistics
################################
################################

# Variables and labels
varz.s <- c("start_year","POW_max","MIA_max","KIA_max","WIA_max","BATTLE_SIZE","size_ratio","polity2_new","gr_polity2_new","cinc","diff_cinc","ler1","ler2","ldr_surrendered","recruit","opp_geneva","DEPLOY_DIST","initiator")
labz.s <- c("Start Year","POW (personnel)","MIA (personnel)","KIA (personnel)","WIA (personnel)","Battle size","Force ratio","Polity 2","More democratic","CINC","More powerful","LER side A","LER side B","Commander surrendered","Professional army","Geneva (opponent)","Deployment distance","Initiator")

# Make table
sum.tab <- data.frame(
  Variable=labz.s,
  Mean=apply(dataX[,varz.s],2,function(x){mean(x,na.rm=T)}),
  SD=apply(dataX[,varz.s],2,function(x){sd(x,na.rm=T)}),
  Min=apply(dataX[,varz.s],2,function(x){min(x,na.rm=T)}),
  Max=apply(dataX[,varz.s],2,function(x){max(x,na.rm=T)}),
  Missing=apply(dataX[,varz.s],2,function(x){sum(is.na(x),na.rm=T)})
)

# Save
print.xtable(xtable(sum.tab,digits = 2),file = "Output/Tables/tab_1_sumtab.tex",include.rownames = FALSE)


#################################
#################################
## Figure 1: map & timeline
#################################
#################################

data("wrld_simpl")

# plot order
areaz <- c()
for(i in 1:length(locations.poly)){areaz[i] <- locations.poly@polygons[[i]]@Polygons[[1]]@area}
# map
png("Output/Figures/fig_1_map.png",width=8,height=4,units="in",res=300)
par(mar=c(0,0,0,0))
plot(wrld_simpl[wrld_simpl$POP2005>50000,],col="bisque",border="white",lwd=1)
for(i in 1:length(areaz)){plot(locations.poly[rev(order(areaz))[i],],add=T,col=rgb(0,.5,.5,.5),border="black",lwd=.5)}
legend("bottomright",fill=rgb(0,.5,.5,.75),legend="Battle locations, 1939-2011",bty="n",cex=1.5)
dev.off()


# timeline
ticker <- data.frame(YEAR=rep(1939:2015,each=12),MONTH=rep(1:12,length(1939:2015)))
ticker$YRMO <- ticker$YEAR*100+ticker$MONTH
ticker$MID <- 1:nrow(ticker)
bidz <- sort(unique(dataX$BATTLE_ID))
ticker$BATTLES <- 0
i <- 598
for(i in 1:length(bidz)){print(i)
  subz <- dataX[which(dataX$BATTLE_ID%in%bidz[i])[1],]
  subz[,1:20]
  if(!is.na(subz$end_yrmo)){ticker[which(ticker$YRMO%in%subz$start_yrmo):(which(ticker$YRMO%in%subz$end_yrmo)+1),"BATTLES"] <- ticker[which(ticker$YRMO%in%subz$start_yrmo):(which(ticker$YRMO%in%subz$end_yrmo)+1),"BATTLES"]+1}
  if(is.na(subz$end_yrmo)){ticker[which(ticker$YRMO%in%subz$start_yrmo):which(ticker$YRMO%in%subz$start_yrmo)+1,"BATTLES"] <- ticker[which(ticker$YRMO%in%subz$start_yrmo):which(ticker$YRMO%in%subz$start_yrmo)+1,"BATTLES"]+1}
}

# plot
png("Output/Figures/fig_1_time.png",width=8,height=2,units="in",res=300)
par(mar=c(2,4,.5,.5))
plot(0,0,col=NA,xlim=range(ticker$MID),ylim=range(ticker$BATTLES),bty="n",xaxt="n",ylab="Battles per month",cex.axis=.8)
axis(1,at=ticker$MID[seq(1,max(ticker$MID),by=12)],labels=NA)
axis(1,at=ticker$MID[seq(13,max(ticker$MID),by=12*5)],labels=ticker$YEAR[seq(13,max(ticker$MID),by=12*5)],lwd=2)
segments(x0=ticker$MID,x1=ticker$MID,y0=0,y1=ticker$BATTLES,col=rgb(0,.5,.5,1))
dev.off()


#################################
#################################
## Figure 2 (OLD): 3-D Network plot
#################################
#################################

datez <- dataX[,c("BATTLE_ID","start_tid","end_tid","start_date","end_date")]
datez <- datez[!duplicated(datez$BATTLE_ID),]
datez <- datez[order(datez$BATTLE_ID),]
locz <- merge(locations.poly,datez,by.x="BATTLE_ID",by.y="BATTLE_ID",all.x=T,all.y=F)
range0 <- bbox(locz)[2,]

# time matrix
startdate <- min(datez$start_date,na.rm=T)
enddate <- max(datez$end_date,na.rm=T)
startyr <- as.numeric(substr(startdate,1,4))
startmo <- as.numeric(substr(startdate,5,6))
startdy <- as.numeric(substr(startdate,7,8))
ndays <- function(mo,yr){
  out <- ifelse(mo%in%c(1,3,5,7,8,10,12),31,ifelse(mo%in%c(4,6,9,11),30,ifelse(yr%%4==0,29,28)))
  return(out)}
ymax <- as.numeric(max(substr(datez$end_date,1,4),na.rm=T))-as.numeric(min(substr(datez$start_date,1,4),na.rm=T))+1
ticker <- lapply(0:ymax,function(y){
  days <- c()
  for(i in 1:12){days <- c(days,1:ndays(1:12,startyr+y)[i])}
  mos <- rep(1:12,ndays(1:12,startyr+y))
  yrs <- rep(startyr+y,length(days))
  data.frame(YEAR=yrs,MONTH=mos,DAY=days)
})
ticker <- do.call(rbind,ticker)
ticker$DATE <- ticker$YEAR*10000+ticker$MONTH*100+ticker$DAY
ticker <- ticker[ticker$DATE>=startdate&ticker$DATE<=enddate,]
ticker$TID <- 1:nrow(ticker)
head(ticker)

# basemap
data("wrld_simpl")
map <- wrld_simpl[wrld_simpl[["POP2005"]]>50000,]
r <- raster(ncol=100, nrow=100)
map <- rasterize(map, r, fun="last")
map[!is.na(map)] <- 1

# connections
disag <- sort(unique(dataX$COW.number))
x0 <- 1
r <- 0.01
k <- 1
i < 10
links.list <- lapply(1:length(disag),function(k){print(k)
  subz <- dataX[which(dataX$COW.number==disag[k]),]
  subz <- subz[order(subz$start_tid),]
  row.names(subz) <- 1:nrow(subz)
  # samei
  links_samei <- lapply(2:nrow(subz),function(i){#print(i)
    ix <- intersect(which(subz[,"code"]%in%subz[i,"code"]),1:(i-1))
    if(length(ix)>0){
      x <- x0
      td <- subz[i,"start_tid"]-subz[ix,"start_tid"]
      data.frame(war=disag[k],type="samei",code=subz[i,"code"],battle_id1=subz[i,"BATTLE_ID"],battle_id2=subz[ix,"BATTLE_ID"],weight=x/(1+r)^td)
    }
  })
  links_samei <- do.call(rbind,links_samei)
  # samej
  links_samej <- lapply(2:nrow(subz),function(i){#print(i)
    opp.code <- subz[which(subz$BATTLE_ID%in%subz[i,"BATTLE_ID"]&subz$initiator%in%1-subz[i,"initiator"]),"code"]
    opp.bat.id <- subz[subz$code%in%opp.code,"BATTLE_ID"]
    ix <- intersect(which(subz$BATTLE_ID%in%opp.bat.id&(!subz$code%in%opp.code)),1:(i-1))
    if(length(ix)>0){
      x <- x0
      td <- subz[i,"start_tid"]-subz[ix,"start_tid"]
      data.frame(war=disag[k],type="samej",code=subz[i,"code"],battle_id1=subz[i,"BATTLE_ID"],battle_id2=subz[ix,"BATTLE_ID"],weight=x/(1+r)^td)
    }
  })
  links_samej <- do.call(rbind,links_samej)
  links.mat <- rbind(links_samei,links_samej)
  links.mat
})
links.mat <- do.call(rbind,links.list)
# add x,y values
locz.xy <- data.frame(battle_id=locz[["BATTLE_ID"]],long=NA,lat=NA)
for(i in 1:length(locz)){
  coordz <- locz[i,]@polygons[[1]]@Polygons[[1]]@coords
  locz.xy[i,c("long","lat")] <- apply(coordz,2,mean)
}
locz.xy1 <- locz.xy; names(locz.xy1) <- paste0(names(locz.xy1),"1")
locz.xy2 <- locz.xy; names(locz.xy2) <- paste0(names(locz.xy2),"2")
links.mat <- merge(links.mat,locz.xy1,by="battle_id1",all.x=T,all.y=F)
links.mat <- merge(links.mat,locz.xy2,by="battle_id2",all.x=T,all.y=F)
# add z values
datez1 <- datez[,c("BATTLE_ID","start_tid")]; names(datez1) <- tolower(paste0(names(datez1),"1"))
datez2 <- datez[,c("BATTLE_ID","start_tid")]; names(datez2) <- tolower(paste0(names(datez2),"2"))
links.mat <- merge(links.mat,datez1,by="battle_id1",all.x=T,all.y=F)
links.mat <- merge(links.mat,datez2,by="battle_id2",all.x=T,all.y=F)
head(links.mat)
links.mat <- links.mat[!is.na(links.mat$weight),]


##
# plot (WWII battles)
##
# parameters
i <- 10
theta0 <- 40
phi0 <- 15
r0 <- 3
d0 <- 1
maxt <- ticker[ticker$DATE==19460101,"TID"]
zlim0 <- c(0,maxt)
res <- persp(map,zlim=zlim0,plot=F,theta=theta0,phi=phi0,r=r0,d=d0,ylim=range0,scale=T)
type0 <- "samej"

for(o in 1:2){
  type0 <- c("samei","samej")[o]
  png(file=paste0("Output/Figures/fig_2X_net_",type0,"_wwii.png"),width=8,height=8,res=300,units="in")
  par(mar=c(0,0,0,0))
  persp(map,zlim=zlim0,col="bisque",border=NA,theta=theta0,phi=phi0,r=r0,d=d0,zlab="Time",xlab="",ylab="",ylim=range0,axes=F,box=T,lwd=.5,scale=T)
  # z axis
  min.x <- bbox(locz)[1,1]
  min.y <- bbox(locz)[2,1]
  z.axis <- ticker[which(ticker$YEAR%%1==0&ticker$MONTH==1&ticker$DAY==1&ticker$DATE<19460101),"TID"]
  tick.start <- trans3d(min.x, min.y, z.axis, res)
  tick.end <- trans3d(min.x, (min.y - 1), z.axis, res)
  segments(tick.start$x, tick.start$y, tick.end$x, tick.end$y)
  z.labels <- ticker[which(ticker$YEAR%%1==0&ticker$MONTH==1&ticker$DAY==1&ticker$DATE<19460101),"YEAR"]
  label.pos <- trans3d(min.x, (min.y - 3), z.axis, res)
  text(label.pos$x, label.pos$y, labels=z.labels, adj=c(1, NA), cex=0.7)
  # # draw links
  if(type0=="samei"){
    links.mati <- links.mat[links.mat$type=="samei"&links.mat$start_tid1<maxt,]
    link.start <- trans3d(links.mati$long1, links.mati$lat1, links.mati$start_tid1, res)
    link.end <- trans3d(links.mati$long2, links.mati$lat2, links.mati$start_tid2, res)
    segments(link.start$x, link.start$y, link.end$x, link.end$y,col=rgb(.8,0,.1,alpha=links.mati$weight/2),lwd=.5)
  }
  if(type0!="samei"){
    links.matj <- links.mat[links.mat$type=="samej"&links.mat$start_tid1<maxt,]
    link.start <- trans3d(links.matj$long1, links.matj$lat1, links.matj$start_tid1, res)
    link.end <- trans3d(links.matj$long2, links.matj$lat2, links.matj$start_tid2, res)
    segments(link.start$x, link.start$y, link.end$x, link.end$y,col=rgb(.8,0,.1,alpha=links.matj$weight/2),lwd=.5)
  }
  #battles
  locz. <- locz[as.data.frame(locz)[,"start_tid"]<maxt,]
  for(i in 1:length(locz.)){
    coordz <- locz.[i,]@polygons[[1]]@Polygons[[1]]@coords
    t <- as.data.frame(locz.)[i,"start_tid"]
    ptz <- trans3d(x=coordz[,1],y=coordz[,2],z=t,pmat=res)
    polygon(ptz,col=rgb(0,.5,.5,.5),border="black",lwd=.4)
  }
  dev.off()
}



#################################
#################################
## Tables 2 and 3: main analyses
#################################
#################################

dvarz <- "lPOW_max"

# Model 1 (diffusion + FE)
form1a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samei_d_",tolower(dvarz)))
form1b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samej_d_",tolower(dvarz)))
form1c <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_theateri_d_",tolower(dvarz)))
form1d <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_theaterj_d_",tolower(dvarz)))

# Model 3 (diffusion + covariates + FE + polity NA fix)
form2a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samei_d_",tolower(dvarz)))
form2b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samej_d_",tolower(dvarz)))
form2c <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_theateri_d_",tolower(dvarz)))
form2d <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_theaterj_d_",tolower(dvarz)))

# size ratio
form3a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+lsize_ratio_2+st_w_samei_d_",tolower(dvarz)))
form3b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+lsize_ratio_2+st_w_samej_d_",tolower(dvarz)))
form3c <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+lsize_ratio_2+st_w_theateri_d_",tolower(dvarz)))
form3d <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+lsize_ratio_2+st_w_theaterj_d_",tolower(dvarz)))

## Run models
x0 <- 1e-4  # lower bound for r
x1 <- .1   # upper bound for r
x0. <- .01  # lower bound for p
x1. <- 1   # upper bound for p
runs <- 50  # number of iterations

# find optimal discount rate via AIC
formz <- form1a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1a <- mod
formz <- form1b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1b <- mod
formz <- form1c
pbest <- optimDiscount_geo(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly,runs=runs,t_bat=FALSE); pbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1c <- mod
formz <- form1d
pbest <- optimDiscount_geo(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly,runs=runs,t_bat=FALSE); pbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1d <- mod
formz <- form2a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2a <- mod
formz <- form2b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2b <- mod
formz <- form2c
pbest <- optimDiscount_geo(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly,runs=runs,t_bat=FALSE); pbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2c <- mod
formz <- form2d
summary(dataX.W$battle_size)
pbest <- optimDiscount_geo(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly,runs=runs,t_bat=FALSE); pbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2d <- mod
formz <- form3a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod3a <- mod
formz <- form3b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod3b <- mod
formz <- form3c
pbest <- optimDiscount_geo(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly,runs=runs,t_bat=FALSE); pbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod3c <- mod
formz <- form3d
summary(dataX.W$battle_size)
pbest <- optimDiscount_geo(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly,runs=runs,t_bat=FALSE); pbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod3d <- mod


# print table (time)
modz <- list(mod1a,mod1b,mod2a,mod2b,mod3a,mod3b)
keepz <- c("st_ldeploy_dist","initiator","st_diff_cinc","recruit","gr_polity2_new","opp_geneva","start_year","st_battle_size","lsize_ratio_2","st_w_samei_d_lpow_max","st_w_samej_d_lpow_max")
labz <- c("Deployment distance","Initiator","More powerful","Professional army","More democratic","Geneva (opponent)","Start year","Battle Size","log(force ratio)","W(same combatant)","W(same opponent)")
stargazer(modz,keep=keepz,covariate.labels=labz ,dep.var.labels = "log(POWs)",ci=T,out="Output/Tables/tab_2_main_time.tex",keep.stat=c("n","ll","aic","ubre","adj.rsq"),star.cutoffs = c(.05,.01,.001),digits = 2,add.lines = list(c("RMSE",round(c(rmse(mod1a),rmse(mod1b),rmse(mod2a),rmse(mod2b),rmse(mod3a),rmse(mod3b)),2)),c("AIC",round(c(mod1a$aic,mod1b$aic,mod2a$aic,mod2b$aic,mod3a$aic,mod3b$aic),2))))


# print table (space)
modz <- list(mod1c,mod1d,mod2c,mod2d,mod3c,mod3d)
keepz <- c("st_ldeploy_dist","initiator","st_diff_cinc","recruit","gr_polity2_new","opp_geneva","start_year","st_battle_size","lsize_ratio_2","st_w_theateri_d_lpow_max","st_w_theaterj_d_lpow_max")
labz <- c("Deployment distance","Initiator","More powerful","Professional army","More democratic","Geneva (opponent)","Start year","Battle Size","log(force ratio)","W(same combatant)","W(same opponent)")
stargazer(modz,keep=keepz,covariate.labels=labz ,dep.var.labels = "log(POWs)",ci=T,out="Output/Tables/tab_3_main_space.tex",keep.stat=c("n","ll","aic","ubre","adj.rsq"),star.cutoffs = c(.05,.01,.001),digits = 2,add.lines = list(c("RMSE",round(c(rmse(mod1c),rmse(mod1d),rmse(mod2c),rmse(mod2d),rmse(mod3c),rmse(mod3d)),2)),c("AIC",round(c(mod1c$aic,mod1d$aic,mod2c$aic,mod2d$aic,mod3c$aic,mod3d$aic),2))))



##################################
##################################
##################################
## Figure 2: simulations / main results
##################################
##################################
##################################


dvarz <- "POW_max"
form2a <- as.formula(paste0("pow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+w_samei_d_",tolower(dvarz)))
form2b <- as.formula(paste0("pow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+w_samej_d_",tolower(dvarz)))

## Run models
x0 <- 1e-4  # lower bound for r
x1 <- .1    # upper bound for r
runs <- 50  # number of iterations
# find optimal discount rate via AIC
formz <- form2a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2a <- mod
xvar <- "w_samei_d_pow_max"
xs1 <- c(quantile(dataX.W[,xvar],.01,na.rm=T),quantile(dataX.W[,xvar],.99,na.rm=T))
xs1 <- c(0,300000)
predz1 <- predmanGAM(mod=mod2a,var=xvar,vals=xs1,quadratic=F,fax=list(cow.number=139,code=365))
formz <- form2b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2b <- mod
xvar <- "w_samej_d_pow_max"
xs2 <- c(quantile(dataX.W[,xvar],.01,na.rm=T),quantile(dataX.W[,xvar],.99,na.rm=T))
xs1 <- c(0,300000)
predz2 <- predmanGAM(mod=mod2b,var=xvar,vals=xs2,quadratic=F,fax=list(cow.number=139,code=365))

# QOIs
predz1$y
predz1$y[2]-predz1$y[1]
predz1$rr
xs1
predz2$y
predz2$y[2]-predz2$y[1]
predz2$rr
xs2

# Simulations
mod <- mod2a
xvar <- "w_samei_d_pow_max"
xs <- seq(min(dataX.W[,xvar],na.rm=T),max(dataX.W[,xvar],na.rm=T),length.out = 100)
predz <- predmanGAM(mod=mod,var=xvar,vals=xs,quadratic=F,fax=list(cow.number=139,code=365))
ylimz <- range(c(predz$lo,predz$up))
predz0 <- predz
ylimz <- range(c(predz$lo,predz$up))
xs <- seq(min(dataX.W[,xvar],na.rm=T),max(dataX.W[,xvar],na.rm=T),length.out = 100)
mod <- mod2b
xvar <- "w_samej_d_pow_max"
predz <- predmanGAM(mod=mod,var=xvar,vals=xs,quadratic=F,fax=list(cow.number=139,code=365))
ylabz <- "Surrender in current battle (# troops)"
xlabz <- "Surrender in previous battles (# troops)"


# Plot
png("Output/Figures/fig_2_pred_pow.png",width=6,height=4,res=300,units="in")
par(xpd=F,mar=c(4,4,.5,.5))
plot(x=0,y=0,xlim=range(xs),ylim=ylimz,bty="n",xaxt="n",yaxt="n",xlab=xlabz,ylab=ylabz,col=NA)
# rect(xleft=min(xs),xright = max(xs),ybottom = ylimz[1]-abs(ylimz[2]-ylimz[1])/2,ytop = ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",border=NA)
segments(x0=min(xs),x1=max(xs),y0=pretty(c(predz0$lo,predz0$up)),y1=pretty(c(predz0$lo,predz0$up)),col="gray90")
segments(x0=min(xs),x1=max(xs),y0=pretty(c(predz0$lo,predz0$up),n=10),y1=pretty(c(predz0$lo,predz0$up),n=10),col="gray90",lwd=.5)
segments(x0=pretty(min(xs):max(xs)),x1=pretty(min(xs):max(xs)),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90")
segments(x0=pretty(min(xs):max(xs),n=10),x1=pretty(min(xs):max(xs),n=10),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",lwd=.5)
polygon(x=c(xs,rev(xs)),y=c(predz$lo,rev(predz$up)),border=NA,col=rgb(0,0,.7,.25))
lines(x=xs,y=predz$y,lwd=2)
polygon(x=c(xs,rev(xs)),y=c(predz0$lo,rev(predz0$up)),border=NA,col=rgb(.7,0,0,.25))
lines(x=xs,y=predz0$y,lwd=2)
par(xpd=F)
axis(1,at=pretty(xs),labels = gdata:::trim(format(pretty(xs),scientific=F,big.mark=',')),cex.axis=.8,lwd=NA,lwd.ticks = 1)
axis(2,at=pretty(c(predz0$lo,predz0$up)),labels=format(pretty(c(predz0$lo,predz0$up)),scientific=F,big.mark=','),las=1,cex.axis=.7,mgp=c(0,-.25,-1),lwd=NA,lwd.ticks = 1)
legend(x=0,y=ylimz[2]-(ylimz[2]-ylimz[1])/50,bty="n",fill=c(rgb(.7,0,0,.25),rgb(0,0,.7,.25)),legend=c("same army","other armies fighting same opponent"),title="Previous surrender within",title.adj = .1)
dev.off()




#################################
#################################
## Table 4: Geneva interactions
#################################
#################################


dvarz <- "lPOW_max"

# Model 1 (Geneva interaction + FE)
form1a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+opp_geneva+st_w_samei_d_",tolower(dvarz),"+I(opp_geneva*st_w_samei_d_",tolower(dvarz),")"))
form1b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+opp_geneva+st_w_samej_d_",tolower(dvarz),"+I(opp_geneva*st_w_samej_d_",tolower(dvarz),")"))

# Model 3 (Geneva interaction + FE + covariates + polity)
form2a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samei_d_",tolower(dvarz),"+I(opp_geneva*st_w_samei_d_",tolower(dvarz),")"))
form2b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samej_d_",tolower(dvarz),"+I(opp_geneva*st_w_samej_d_",tolower(dvarz),")"))

## Run models
x0 <- 1e-4  # lower bound for r
x1 <- .1   # upper bound for r
x0. <- .01  # lower bound for p
x1. <- 1   # upper bound for p
runs <- 50  # number of iteration
# find optimal discount rate via AIC
formz <- form1a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1a <- mod
formz <- form1b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1b <- mod
formz <- form2a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2a <- mod
formz <- form2b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2b <- mod


# print table (time)
modz <- list(mod1a,mod1b,mod2a,mod2b)
keepz <- c("opp_geneva","st_w_samei_d_lpow_max","I(opp_geneva * st_w_samei_d_lpow_max)","st_w_samej_d_lpow_max","I(opp_geneva * st_w_samej_d_lpow_max)")
labz <- c("Geneva (opponent)","W(same combatant)","Interaction (same combatant)","W(same opponent)","Interaction (same opponent)")
stargazer(modz,keep=keepz,covariate.labels=labz ,dep.var.labels = "log(POWs)",ci=T,keep.stat=c("n","ll","aic","ubre","adj.rsq"),star.cutoffs = c(.05,.01,.001),digits = 2,out="Output/Tables/tab_4_geneva.tex",add.lines = list(c("RMSE",round(c(rmse(mod1a),rmse(mod1b),rmse(mod2a),rmse(mod2b)),2)),c("AIC",round(c(mod1a$aic,mod1b$aic,mod2a$aic,mod2b$aic),2))))



#################################
#################################
## Table 5: surrender by commanders
#################################
#################################

# Model 1 (troops ~ commanders + FE)
dvarz <- "ldr_surrendered"
form1a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samei_d_",tolower(dvarz)))
form1b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samej_d_",tolower(dvarz)))

# Model 1 (troops ~ commanders + FE + covariates)
dvarz <- "ldr_surrendered"
form2a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samei_d_",tolower(dvarz)))
form2b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samej_d_",tolower(dvarz)))

## Run models
x0 <- 1e-4  # lower bound for r
x1 <- .1   # upper bound for r
x0. <- 1e-2  # lower bound for p
x1. <- 1   # upper bound for p
runs <- 50  # number of iterations
# find optimal discount rate via AIC
formz <- form1a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1a <- mod
formz <- form1b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1b <- mod
formz <- form2a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2a <- mod
formz <- form2b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2b <- mod


# Model 3 (commanders ~ commanders + FE)
dvarz <- "ldr_surrendered"
form3a <- as.formula(paste0("ldr_surrendered~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samei_d_",tolower(dvarz)))
form3b <- as.formula(paste0("ldr_surrendered~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samej_d_",tolower(dvarz)))

# Model 4 (commanders ~ commanders + FE + covariates)
dvarz <- "ldr_surrendered"
form4a <- as.formula(paste0("ldr_surrendered~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+opp_geneva+start_year+st_battle_size+st_w_samei_d_",tolower(dvarz)))
form4b <- as.formula(paste0("ldr_surrendered~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+opp_geneva+start_year+st_battle_size+st_w_samej_d_",tolower(dvarz)))

## Run models
x0 <- 1e-4  # lower bound for r
x1 <- .1   # upper bound for r
x0. <- 1e-2  # lower bound for p
x1. <- 1   # upper bound for p
runs <- 50  # number of iterations
# find optimal discount rate via AIC
formz <- form3a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="binomial"); summary(mod); mod3a <- mod
formz <- form3b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="binomial"); summary(mod); mod3b <- mod
formz <- form4a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="binomial"); summary(mod); mod4a <- mod
formz <- form4b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="binomial"); summary(mod); mod4b <- mod

# print table
modz <- list(mod1a,mod1b,mod2a,mod2b,mod3a,mod3b,mod4a,mod4b)
keepz <- c("st_w_samei_d_ldr_surrendered","st_w_samej_d_ldr_surrendered","st_w_samei_d_lpow_max")
labz <- c("W(same combatant)","W(same opponent)")
stargazer(modz,keep=keepz,covariate.labels=labz,dep.var.labels = c("log(POWs)","Commander"),ci=T,out="Output/Tables/tab_5_pa.tex",keep.stat=c("n","ll","aic","ubre","adj.rsq"),star.cutoffs = c(.05,.01,.001),digits = 2)
statz <- rbind(c(rmse(mod1a),rmse(mod1b),rmse(mod2a),rmse(mod2b),rmse(mod3a),rmse(mod3b),rmse(mod4a),rmse(mod4b)),c(mod1a$aic,mod1b$aic,mod2a$aic,mod2b$aic,mod3a$aic,mod3b$aic,mod4a$aic,mod4b$aic))
row.names(statz) <- c("RMSE","AIC")
print.xtable(xtable(statz),digits=3,include.colnames = F,file="Output/Tables/tab_5_pa_statz.tex")

# y-standardization
st.coef <- matrix("",length(keepz),length(modz)+1);
st.coef[,1] <- labz
st.ci <- matrix("",length(keepz),length(modz)+1);st.ci[,1] <- ""
num <- 3
k <- 4
for(k in 1:length(modz)){
  if(modz[[k]]$family[1]!="binomial"){
    st.b <- summary(modz[[k]])$p.table[,1]
    starzy <- ifelse(summary(modz[[k]])$p.table[,4]<.001,"***",ifelse(summary(modz[[k]])$p.table[,4]<.01,"**",ifelse(summary(modz[[k]])$p.table[,4]<.05,"*","")))
    st.lo <- (summary(modz[[k]])$p.table[,1]-1.96*summary(modz[[k]])$p.table[,2])
    st.up <- (summary(modz[[k]])$p.table[,1]+1.96*summary(modz[[k]])$p.table[,2])
    st.coef[keepz%in%names(st.b),k+1] <- paste0(round(st.b[keepz[keepz%in%names(st.b)]],num),starzy[keepz[keepz%in%names(st.b)]])
    st.ci[keepz%in%names(st.b),k+1] <- paste0("(",round(st.lo,num),", ",round(st.up,num),")")[names(st.b)%in%keepz[keepz%in%names(st.b)]]
  }
  if(modz[[k]]$family[1]=="binomial"){
    st.b <- summary(modz[[k]])$p.table[,1]/sd(modz[[k]]$y)
    starzy <- ifelse(summary(modz[[k]])$p.table[,4]<.001,"***",ifelse(summary(modz[[k]])$p.table[,4]<.01,"**",ifelse(summary(modz[[k]])$p.table[,4]<.05,"*","")))
    st.lo <- (summary(modz[[k]])$p.table[,1]-1.96*summary(modz[[k]])$p.table[,2])/sd(modz[[k]]$y)
    st.up <- (summary(modz[[k]])$p.table[,1]+1.96*summary(modz[[k]])$p.table[,2])/sd(modz[[k]]$y)
    st.coef[keepz%in%names(st.b),k+1] <- paste0(round(st.b[keepz[keepz%in%names(st.b)]],num),starzy[keepz[keepz%in%names(st.b)]])
    st.ci[keepz%in%names(st.b),k+1] <- paste0("(",round(st.lo,num),", ",round(st.up,num),")")[names(st.b)%in%keepz[keepz%in%names(st.b)]]
  }
}
st.table <- matrix("",length(keepz)*3,length(modz)+1)
st.table[which(1:nrow(st.table)%%3==0)-2,] <- st.coef
st.table[which(1:nrow(st.table)%%3==0)-1,] <- st.ci
print.xtable(xtable(st.table),digits=3,include.colnames = F,file="Output/Tables/tab_5_pa_stcoef.tex",include.rownames=F)






#################################
#################################
## Table 6: placebo tests / military effectiveness
#################################
#################################

# Model 1 (KIA+WIA + FE + covariates)
dvarz <- "ltotal_cas1"
form1a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samei_d_",tolower(dvarz)))
form1b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samej_d_",tolower(dvarz)))

## Run models
x0 <- 1e-4  # lower bound for r
x1 <- .1   # upper bound for r
runs <- 50  # number of iterations
# find optimal discount rate via AIC
formz <- form1a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1a <- mod
formz <- form1b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod1b <- mod


# Model 2 (LER + FE + covariates)
dvarz <- "ler1"
form2a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samei_d_",tolower(dvarz)))
form2b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samej_d_",tolower(dvarz)))

## Run models
x0 <- 1e-4  # lower bound for r
x1 <- .1   # upper bound for r
runs <- 50  # number of iterations

# find optimal discount rate via AIC
formz <- form2a
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2a <- mod
formz <- form2b
rbest <- optimDiscount(formz=formz,X=dataX,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod); mod2b <- mod

# print table
modz <- list(mod1a,mod1b,mod2a,mod2b)
keepz <- c("st_w_samei_d_ltotal_cas1","st_w_samej_d_ltotal_cas1","st_w_samei_d_ler1","st_w_samej_d_ler1")
labz <- c("W(same combatant): KIA+WIA","W(same opponent): KIA+WIA","W(same combatant): LER","W(same opponent): LER")
stargazer(modz,keep=keepz,covariate.labels=labz,dep.var.labels = "log(POWs)",ci=T,out="Output/Tables/tab_6_placebo.tex",keep.stat=c("n","ll","aic","ubre","adj.rsq"),star.cutoffs = c(.05,.01,.001),digits = 2)
statz <- rbind(c(rmse(mod1a),rmse(mod1b),rmse(mod2a),rmse(mod2b)),c(mod1a$aic,mod1b$aic,mod2a$aic,mod2b$aic))
row.names(statz) <- c("RMSE","AIC")
print.xtable(xtable(statz),digits=3,include.colnames = F,file="Output/Tables/tab_6_placebo_statz.tex")






#################################
#################################
## Figure 3: Sensitivity analysis 1: discount rate
#################################
#################################

dvarz <- "lPOW_max"
# Model 1 (diffusion + FE)
form1a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samei_d_",tolower(dvarz)))
form1b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samej_d_",tolower(dvarz)))
form1c <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_theateri_d_",tolower(dvarz)))
form1d <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_theaterj_d_",tolower(dvarz)))
rz <- 1:1000/1000
r <- 2
rz[r]
write.file <- 'R_progress.txt'

# Loop over rz vector (NOTE: this takes a long time to run)
rlist <- mclapply(seq_along(rz),function(r){print(rz[r])

  # monitor progress
  file.create(write.file)
  fileConn<-file(write.file)
  writeLines(paste0(r,'/',length(rz),' ',rz[r]), fileConn)
  close(fileConn)
  # monitor from a console with
  # tail -c +0 -f ~/Dropbox/BATTLEDATA/R_progress.txt

  # Pick r
  rbest <- rz[r]

  # Fit models
  formz <- form1a
  dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
  mod <- glm(formz,data=dataX.W,family="gaussian");# summary(mod)
  names(summary(mod))
  summary(mod)$coefficients
  beta_1 <- summary(mod)$coefficients[paste0("st_w_samei_d_",tolower(dvarz)),1]
  se_1 <- summary(mod)$coefficients[paste0("st_w_samei_d_",tolower(dvarz)),2]
  aic_1 <- mod$aic
  formz <- form1b
  dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
  mod <- glm(formz,data=dataX.W,family="gaussian");# summary(mod)
  beta_2 <- summary(mod)$coefficients[paste0("st_w_samej_d_",tolower(dvarz)),1]
  se_2 <- summary(mod)$coefficients[paste0("st_w_samej_d_",tolower(dvarz)),2]
  aic_2 <- mod$aic
  formz <- form1c
  dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,p=rbest,polyz=locations.poly,t_bat=FALSE)
  mod <- glm(formz,data=dataX.W,family="gaussian");# summary(mod)
  beta_3 <- summary(mod)$coefficients[paste0("st_w_theateri_d_",tolower(dvarz)),1]
  se_3 <- summary(mod)$coefficients[paste0("st_w_theateri_d_",tolower(dvarz)),2]
  aic_3 <- mod$aic
  formz <- form1d
  dataX.W <- varyDiscount2(X=dataX,dvarz=dvarz,xvarz=xvarz1,p=rbest,polyz=locations.poly,t_bat=FALSE)
  mod <- glm(formz,data=dataX.W,family="gaussian")#; summary(mod)
  beta_4 <- summary(mod)$coefficients[paste0("st_w_theaterj_d_",tolower(dvarz)),1]
  se_4 <- summary(mod)$coefficients[paste0("st_w_theaterj_d_",tolower(dvarz)),2]
  aic_4 <- mod$aic

  # Extract QOIs
  rmat <- data.frame(r=rbest,beta_1=beta_1,se_1=se_1,aic_1=aic_1,beta_2=beta_2,se_2=se_2,aic_2=aic_2,beta_3=beta_3,se_3=se_3,aic_3=aic_3,beta_4=beta_4,se_4=se_4,aic_4=aic_4)
  rmat
},mc.preschedule = FALSE, mc.set.seed = TRUE,mc.silent = FALSE, mc.cores = detectCores())
rmatz <- do.call(rbind,rlist)
# save(rmatz,file="Output/RobustCheck_DiscountRate.RData")
load("Output/RobustCheck_DiscountRate.RData")


################
## Plot
################

dev.off()

png(file="Output/Figures/fig_3_robust_discount.png",height=5,width=8.5,res=300,units="in")
par(mar=c(4,4,0.5,0.5),xpd=0,mfrow=c(2,2))

# Plot (Model 1)
ylimz <- range(pretty(c(-.05,rmatz[,"beta_1"]+2*rmatz[,"se_1"],rmatz[,"beta_1"]-2*rmatz[,"se_1"])))
plot(0,0,xlim=c(0,1),ylim=ylimz,col=NA,bty="n",xlab="Discount rate (r), temporal weights",ylab=expression(paste("Effect of past surrender ( ", hat(rho), " )")),xaxt="n",yaxt="n",mgp=c(2.5,1,0))
# rect(xleft=0,xright = 1,ybottom = ylimz[1]-abs(ylimz[2]-ylimz[1])/2,ytop = ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",border=NA)
segments(x0=min(rz),x1=max(rz),y0=pretty(ylimz),y1=pretty(ylimz),col="gray90")
segments(x0=min(rz),x1=max(rz),y0=pretty(ylimz,n=10),y1=pretty(ylimz,n=10),col="gray90",lwd=.5)
segments(x0=pretty(rz),x1=pretty(rz),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90")
segments(x0=pretty(rz,n=10),x1=pretty(rz,n=10),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",lwd=.5)
segments(x0=0,x1=1,y0=0,y1=0,col="gray90",lwd=5)
polygon(x=c(rz,rev(rz)),y=c(rmatz$beta_1-1.96*rmatz$se_1,rev(rmatz$beta_1+1.96*rmatz$se_1)),border=NA,col=rgb(.7,0,0,.25))
lines(x=rz,y=rmatz$beta_1,lwd=2)
axis(1,at=pretty(rz),labels = gdata:::trim(format(pretty(rz),scientific=F,big.mark=',')),cex.axis=.7,lwd=NA,lwd.ticks = 1)
axis(2,at=pretty(ylimz),labels=format(pretty(ylimz),scientific=F,big.mark=','),las=1,cex.axis=.7,mgp=c(3,1,-0.75),lwd=NA,lwd.ticks = 1)
r0 <- rmatz[rmatz$aic_1==min(rmatz$aic_1),"r"]
aic0 <- rmatz[rmatz$aic_1==min(rmatz$aic_1),"aic_1"]
segments(x0=r0,x1=r0,y0=ylimz[1],y1=ylimz[2]-(ylimz[2]-ylimz[1])/10,lty=1,col=rgb(.3,.4,.3,1))
text(x=r0,y=ylimz[1]+(ylimz[2]-ylimz[1])/20,labels=paste0("r*=",r0),pos=4,col=rgb(.3,.4,.3,1))
text(x=1,y=ylimz[2]-(ylimz[2]-ylimz[1])/20,labels=expression(Past~surrender~within~bold(same)~bold(army)),pos=2,cex=1)


# Plot (Model 2)
ylimz <- range(pretty(c(-.05,rmatz[,"beta_2"]+2*rmatz[,"se_2"],rmatz[,"beta_2"]-2*rmatz[,"se_2"])))
plot(0,0,xlim=c(0,1),ylim=ylimz,col=NA,bty="n",xlab="Discount rate (r), temporal weights",ylab=expression(paste("Effect of past surrender ( ", hat(rho), " )")),xaxt="n",yaxt="n",mgp=c(2.5,1,0))
# rect(xleft=0,xright = 1,ybottom = ylimz[1]-abs(ylimz[2]-ylimz[1])/2,ytop = ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",border=NA)
segments(x0=min(rz),x1=max(rz),y0=pretty(ylimz),y1=pretty(ylimz),col="gray90")
segments(x0=min(rz),x1=max(rz),y0=pretty(ylimz,n=10),y1=pretty(ylimz,n=10),col="gray90",lwd=.5)
segments(x0=pretty(rz),x1=pretty(rz),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90")
segments(x0=pretty(rz,n=10),x1=pretty(rz,n=10),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",lwd=.5)
segments(x0=0,x1=1,y0=0,y1=0,col="gray90",lwd=5)
polygon(x=c(rz,rev(rz)),y=c(rmatz$beta_2-1.96*rmatz$se_2,rev(rmatz$beta_2+1.96*rmatz$se_2)),border=NA,col=rgb(0,0,.7,.25))
lines(x=rz,y=rmatz$beta_2,lwd=2)
axis(1,at=pretty(rz),labels = gdata:::trim(format(pretty(rz),scientific=F,big.mark=',')),cex.axis=.7,lwd=NA,lwd.ticks = 1)
axis(2,at=pretty(ylimz),labels=format(pretty(ylimz),scientific=F,big.mark=','),las=1,cex.axis=.7,mgp=c(3,1,-0.75),lwd=NA,lwd.ticks = 1)
r0 <- rmatz[rmatz$aic_2==min(rmatz$aic_2),"r"]
aic0 <- rmatz[rmatz$aic_2==min(rmatz$aic_2),"aic_2"]
segments(x0=r0,x1=r0,y0=ylimz[1],y1=ylimz[2]-(ylimz[2]-ylimz[1])/10,lty=1,col=rgb(.3,.4,.3,1))
text(x=r0,y=ylimz[1]+(ylimz[2]-ylimz[1])/20,labels=paste0("r*=",r0),pos=4,col=rgb(.3,.4,.3,1))
text(x=1,y=ylimz[2]-(ylimz[2]-ylimz[1])/20,labels=expression("..."~within~bold(other)~bold(armies)~bold(fighting)~bold(same)~bold(opponent)),pos=2,cex=1)


# Plot (Model 3)
ylimz <- range(pretty(c(-.05,rmatz[,"beta_3"]+2*rmatz[,"se_3"],rmatz[,"beta_3"]-2*rmatz[,"se_3"])))
plot(0,0,xlim=c(0,1),ylim=ylimz,col=NA,bty="n",xlab="Discount rate (r), geographic weights",ylab=expression(paste("Effect of past surrender ( ", hat(rho), " )")),xaxt="n",yaxt="n",mgp=c(2.5,1,0))
# rect(xleft=0,xright = 1,ybottom = ylimz[1]-abs(ylimz[2]-ylimz[1])/2,ytop = ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",border=NA)
segments(x0=min(rz),x1=max(rz),y0=pretty(ylimz),y1=pretty(ylimz),col="gray90")
segments(x0=min(rz),x1=max(rz),y0=pretty(ylimz,n=10),y1=pretty(ylimz,n=10),col="gray90",lwd=.5)
segments(x0=pretty(rz),x1=pretty(rz),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90")
segments(x0=pretty(rz,n=10),x1=pretty(rz,n=10),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",lwd=.5)
segments(x0=0,x1=1,y0=0,y1=0,col="gray90",lwd=5)
polygon(x=c(rz,rev(rz)),y=c(rmatz$beta_3-1.96*rmatz$se_3,rev(rmatz$beta_3+1.96*rmatz$se_3)),border=NA,col=rgb(.7,0,0,.25))
lines(x=rz,y=rmatz$beta_3,lwd=2)
axis(1,at=pretty(rz),labels = gdata:::trim(format(pretty(rz),scientific=F,big.mark=',')),cex.axis=.7,lwd=NA,lwd.ticks = 1)
axis(2,at=pretty(ylimz),labels=format(pretty(ylimz),scientific=F,big.mark=','),las=1,cex.axis=.7,mgp=c(3,1,-0.75),lwd=NA,lwd.ticks = 1)
r0 <- rmatz[rmatz$aic_3==min(rmatz$aic_3),"r"]
aic0 <- rmatz[rmatz$aic_3==min(rmatz$aic_3),"aic_3"]
segments(x0=r0,x1=r0,y0=ylimz[1],y1=ylimz[2]-(ylimz[2]-ylimz[1])/10,lty=1,col=rgb(.3,.4,.3,1))
text(x=r0,y=ylimz[1]+(ylimz[2]-ylimz[1])/20,labels=paste0("r*=",r0),pos=4,col=rgb(.3,.4,.3,1))
text(x=1,y=ylimz[2]-(ylimz[2]-ylimz[1])/20,labels=expression(Past~surrender~within~bold(same)~bold(army)),pos=2,cex=1)


# Plot (Model 4)
ylimz <- range(pretty(c(-.05,rmatz[,"beta_4"]+2*rmatz[,"se_4"],rmatz[,"beta_4"]-2*rmatz[,"se_4"])))
plot(0,0,xlim=c(0,1),ylim=ylimz,col=NA,bty="n",xlab="Discount rate (r), geographic weights",ylab=expression(paste("Effect of past surrender ( ", hat(rho), " )")),xaxt="n",yaxt="n",mgp=c(2.5,1,0))
# rect(xleft=0,xright = 1,ybottom = ylimz[1]-abs(ylimz[2]-ylimz[1])/2,ytop = ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",border=NA)
segments(x0=min(rz),x1=max(rz),y0=pretty(ylimz),y1=pretty(ylimz),col="gray90")
segments(x0=min(rz),x1=max(rz),y0=pretty(ylimz,n=10),y1=pretty(ylimz,n=10),col="gray90",lwd=.5)
segments(x0=pretty(rz),x1=pretty(rz),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90")
segments(x0=pretty(rz,n=10),x1=pretty(rz,n=10),y0=ylimz[1]-abs(ylimz[2]-ylimz[1])/2,y1=ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",lwd=.5)
segments(x0=0,x1=1,y0=0,y1=0,col="gray90",lwd=5)
polygon(x=c(rz,rev(rz)),y=c(rmatz$beta_4-1.96*rmatz$se_4,rev(rmatz$beta_4+1.96*rmatz$se_4)),border=NA,col=rgb(0,0,.7,.25))
lines(x=rz,y=rmatz$beta_4,lwd=2)
axis(1,at=pretty(rz),labels = gdata:::trim(format(pretty(rz),scientific=F,big.mark=',')),cex.axis=.7,lwd=NA,lwd.ticks = 1)
axis(2,at=pretty(ylimz),labels=format(pretty(ylimz),scientific=F,big.mark=','),las=1,cex.axis=.7,mgp=c(3,1,-0.75),lwd=NA,lwd.ticks = 1)
r0 <- rmatz[rmatz$aic_4==min(rmatz$aic_4),"r"]
aic0 <- rmatz[rmatz$aic_4==min(rmatz$aic_4),"aic_4"]
segments(x0=r0,x1=r0,y0=ylimz[1],y1=ylimz[2]-(ylimz[2]-ylimz[1])/10,lty=1,col=rgb(.3,.4,.3,1))
text(x=r0,y=ylimz[1]+(ylimz[2]-ylimz[1])/20,labels=paste0("r*=",r0),pos=4,col=rgb(.3,.4,.3,1))
text(x=1,y=ylimz[2]-(ylimz[2]-ylimz[1])/20,labels=expression("..."~within~bold(other)~bold(armies)~bold(fighting)~bold(same)~bold(opponent)),pos=2,cex=1)

dev.off()



#################################
#################################
## Figure 4: Japan LER & POW
#################################
#################################


###########################
## LER over time
###########################

head(dataX)
bat.list <- sort(unique(dataX$COW.number))
b <- 1
bat.list[b]
# varz <- c("lPOW_max","lKIA_max","lWIA_max","lMIA_max","ltotal_cas1","ltotal_cas2","ller1","ller2")
# labz <- c("log( POW )","log( KIA )","log( WIA )","log( MIA )","log( Total Casualties )","log( Total Casualties )","log( LER )","log( LER )")
varz <- c("ller2","lPOW_max")
labz <- c("log( LER )","log( POW )")

# Loop over wars
sub.war <- dataX[dataX$COW.number%in%bat.list[b],]

# Time vector
head(sub.war)
date.range <- paste0(substr(range(sub.war$start_yrmo,sub.war$end_yrmo),1,4),"-",substr(range(sub.war$start_yrmo,sub.war$end_yrmo),5,6),"-01")
date.range <- gsub("-|01$","",seq(as.Date(date.range[1]), as.Date(date.range[2]), by="months"))
date.range <- data.frame(TID=seq_along(date.range),YRMO=date.range)
loss.mat <- as.data.frame(matrix(NA,nrow=nrow(date.range),ncol=length(varz)))
names(loss.mat) <- varz
loss.mat <- cbind(date.range,loss.mat)
head(loss.mat)

# Loop over war participants
combtz <- sort(unique(sub.war$code))
k <- 16
combtz.list <- lapply(seq_along(combtz),function(k){
  sub.comb <- sub.war[sub.war$code%in%combtz[k],]
  tail(sub.comb)
  sub.comb$ller2 <- log(sub.comb$ler2+1)
  # Loop over months
  monthz <- loss.mat$YRMO <- as.numeric(as.character(loss.mat$YRMO))
  t <- 40
  for(t in seq_along(monthz)){
    sub.month <- sub.comb[sub.comb$start_yrmo<=monthz[t]&sub.comb$end_yrmo>=monthz[t],]
    if(nrow(sub.month)>0){loss.mat[t,varz] <- colMeans(sub.month[,varz],na.rm=T)}
  }
  list(COW=bat.list[b],CCODE=combtz[k],WAR=sub.comb$conflict[1],COUNTRY=sub.comb$name[1],LOSSES=loss.mat)
})

sapply(combtz.list,function(x){unlist(x["CCODE"])})
which(sapply(combtz.list,function(x){unlist(x["COUNTRY"])})%in%"USSR")

###################
## Plot
###################

# Loop over countries
k <- 1
j <- 1
combtz.n <- c(740)
varz.n <- "ller2"
combtz[k]
dev.off()

for(j in seq_along(varz)){
  png(file=paste0("Output/Figures/fig_4_Japan_",varz[j],".png"),width=6,height=2,res=300,units="in")
  par(mar=c(2,4,.5,.5))
  loss.mat0 <- combtz.list[[which(combtz%in%combtz.n[k])]]$LOSSES
  tail(loss.mat0)
  if(combtz.n[k]%in%c(255,365)){loss.mat0[loss.mat0$YRMO>194505,varz[j]] <- NA}
  ylimz <- range(c(0,loss.mat0[,varz[j]]),na.rm=T)
  par(mar=c(2,4,.5,.5))
  plot(0,0,xlim=c(0,nrow(loss.mat0)),ylim=ylimz,col=NA,bty="n",xlab="",ylab=labz[j],xaxt="n",yaxt="n",mgp=c(2.5,1,0))
  # rect(xleft=0-1,xright = nrow(loss.mat0)+1,ybottom = ylimz[1]-abs(ylimz[2]-ylimz[1])/2,ytop = ylimz[2]+abs(ylimz[2]-ylimz[1])/2,col="gray90",border=NA)
  segments(x0=0,x1=nrow(loss.mat0),y0=pretty(ylimz),y1=pretty(ylimz),col="gray90")
  segments(x0=0,x1=nrow(loss.mat0),y0=pretty(ylimz)-diff(pretty(ylimz))[1]/2,y1=pretty(ylimz)-diff(pretty(ylimz))[1]/2,col="gray90",lwd=.5)
  xs <- 1:nrow(loss.mat0)
  rect(xleft = xs-.4,xright=xs+.4,ybottom = 0,ytop = loss.mat0[,varz[j]],border=NA,col=rgb(.2,.2,.2,.5))
  axis(1,at=loss.mat0$TID[which(substr(loss.mat0$YRMO,5,6)=="01")],labels=substr(loss.mat0$YRMO[which(substr(loss.mat0$YRMO,5,6)=="01")],1,4))
  axis(1,at=loss.mat0$TID,labels=NA,lwd=.5)
  axis(2,at=pretty(ylimz),labels=pretty(ylimz),las=1,mgp=c(3,1,-0.75),lwd=0,lwd.ticks = 1)
  if(varz[j]%in%c("ller1","ller2")){
    segments(x0=-1,x1=nrow(loss.mat0)+1,y0=log(2),y1=log(2),lty=1,col=rgb(1,.2,0,1))
    par(xpd="n")
    if(!combtz.n[k]%in%c(255)){text(x=-3+(nrow(loss.mat0))/30,y=log(2)-.15,labels=paste0("LER = 1:1"),pos=2,col=rgb(1,.2,0,1),cex=.8)}
  }
dev.off()
}



#################################
#################################
## Sensitivity analysis 2: no sea or air battles
#################################
#################################

dvarz <- "lPOW_max"

# Model 1 (diffusion + FE)
form1a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samei_d_",tolower(dvarz)))
form1b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_samej_d_",tolower(dvarz)))
form1c <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_theateri_d_",tolower(dvarz)))
form1d <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_w_theaterj_d_",tolower(dvarz)))

# Model 3 (diffusion + covariates + FE + polity NA fix)
form2a <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samei_d_",tolower(dvarz)))
form2b <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_samej_d_",tolower(dvarz)))
form2c <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_theateri_d_",tolower(dvarz)))
form2d <- as.formula(paste0("st_lpow_max~as.factor(cow.number)+as.factor(code)+spring+summer+fall+st_ldeploy_dist+initiator+st_diff_cinc+recruit+gr_polity2_new+opp_geneva+start_year+st_battle_size+st_w_theaterj_d_",tolower(dvarz)))


## Run models
x0 <- 1e-4  # lower bound for r
x1 <- .1   # upper bound for r
x0. <- .01  # lower bound for p
x1. <- 1   # upper bound for p
runs <- 50  # number of iterations

# find optimal discount rate via AIC
formz <- form1a
rbest <- optimDiscount(formz=formz,X=dataX0,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX0.W <- varyDiscount2(X=dataX0,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX0.W,family="gaussian"); summary(mod); mod1a <- mod
formz <- form1b
rbest <- optimDiscount(formz=formz,X=dataX0,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX0.W <- varyDiscount2(X=dataX0,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX0.W,family="gaussian"); summary(mod); mod1b <- mod
formz <- form1c
pbest <- optimDiscount_geo(formz=formz,X=dataX0,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly0,runs=runs,t_bat=FALSE); pbest
dataX0.W <- varyDiscount2(X=dataX0,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly0,t_bat=FALSE)
mod <- gam(formz,data=dataX0.W,family="gaussian"); summary(mod); mod1c <- mod
formz <- form1d
pbest <- optimDiscount_geo(formz=formz,X=dataX0,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly0,runs=runs,t_bat=FALSE); pbest #
dataX0.W <- varyDiscount2(X=dataX0,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly0,t_bat=FALSE)
mod <- gam(formz,data=dataX0.W,family="gaussian"); summary(mod); mod1d <- mod
formz <- form2a
rbest <- optimDiscount(formz=formz,X=dataX0,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest
dataX0.W <- varyDiscount2(X=dataX0,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX0.W,family="gaussian"); summary(mod); mod2a <- mod
formz <- form2b
rbest <- optimDiscount(formz=formz,X=dataX0,dvarz=dvarz,xvarz=xvarz1,x0=x0,x1=x1,runs=runs,t_bat=FALSE); rbest #
dataX0.W <- varyDiscount2(X=dataX0,dvarz=dvarz,xvarz=xvarz1,r=rbest,t_bat=FALSE)
mod <- gam(formz,data=dataX0.W,family="gaussian"); summary(mod); mod2b <- mod
formz <- form2c
pbest <- optimDiscount_geo(formz=formz,X=dataX0,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly0,runs=runs,t_bat=FALSE); pbest
dataX0.W <- varyDiscount2(X=dataX0,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly0,t_bat=FALSE)
mod <- gam(formz,data=dataX0.W,family="gaussian"); summary(mod); mod2c <- mod
formz <- form2d
pbest <- optimDiscount_geo(formz=formz,X=dataX0,dvarz=dvarz,xvarz=xvarz1,x0=x0.,x1=x1.,polyz=locations.poly0,runs=runs,t_bat=FALSE); pbest #
dataX0.W <- varyDiscount2(X=dataX0,dvarz=dvarz,xvarz=xvarz1,p=pbest,polyz=locations.poly0,t_bat=FALSE)
mod <- gam(formz,data=dataX0.W,family="gaussian"); summary(mod); mod2d <- mod


# print table (time)
modz <- list(mod1a,mod1b,mod2a,mod2b)
keepz <- c("st_ldeploy_dist","initiator","st_diff_cinc","recruit","gr_polity2_new","opp_geneva","start_year","st_battle_size","st_w_samei_d_lpow_max","st_w_samej_d_lpow_max")
labz <- c("Deployment distance","Initiator","More powerful","Professional army","More democratic","Geneva (opponent)","Start year","Battle Size","W(same combatant)","W(same opponent)")
stargazer(modz,keep=keepz,covariate.labels=labz ,dep.var.labels = "log(POWs)",ci=T,out="Output/Tables/tab_2X_nonav_time.tex",keep.stat=c("n","ll","aic","ubre","adj.rsq"),star.cutoffs = c(.05,.01,.001),digits = 2,add.lines = list(c("RMSE",round(c(rmse(mod1a),rmse(mod1b),rmse(mod2a),rmse(mod2b)),2)),c("AIC",round(c(mod1a$aic,mod1b$aic,mod2a$aic,mod2b$aic),2))))


# print table (space)
modz <- list(mod1c,mod1d,mod2c,mod2d)
keepz <- c("st_ldeploy_dist","initiator","st_diff_cinc","recruit","gr_polity2_new","opp_geneva","start_year","st_battle_size","st_w_theateri_d_lpow_max","st_w_theaterj_d_lpow_max")
labz <- c("Deployment distance","Initiator","More powerful","Professional army","More democratic","Geneva (opponent)","Start year","Battle Size","W(same combatant)","W(same opponent)")
stargazer(modz,keep=keepz,covariate.labels=labz ,dep.var.labels = "log(POWs)",ci=T,out="Output/Tables/tab_3X_nonav_space.tex",keep.stat=c("n","ll","aic","ubre","adj.rsq"),star.cutoffs = c(.05,.01,.001),digits = 2,add.lines = list(c("RMSE",round(c(rmse(mod1c),rmse(mod1d),rmse(mod2c),rmse(mod2d)),2)),c("AIC",round(c(mod1c$aic,mod1d$aic,mod2c$aic,mod2d$aic),2))))



